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HI — RECIRCULATING FLOWS 


RESUMO 


Discutem-se os escoamentos bidimensionais em que as 
aproximações de camada limite não são válidas, e os aspec- 
tos matemáticos e físicos do problema. As dificuldades 
intrinsecas do tratamento numérico das equações são exem- 
plificadas, mostrando-se que o aumento do número de Rey- 
nolds conduz a exigências contraditórias no que ae refere a 
precisão e estabilidade. 

Apresentam-se modelos de turbulência e a sua incorpo- 
ração num método de cálculo aplicável a este tepo de escod- 


mentos tal como desenvolvido por D. B. Spalding et ui. 


0 — INTRODUCTION 


In part I [1] the general equations which form 
the actual basis for the study of turbulent 
incompressible flows were deduced and discussed. 
The combination of fundamental equations with 
empiricism or guess, which allowed the proposal 
of turbulence models to be applicable in prediction 
methods of boundary-layer flows, were the object 
of a detailed discussion in part II [2], together 
with the assessment of their practical utility. 

Boundary-layers are the particular class of 
two-dimensional flows where a dominant flow 
direction exists such that the momentum equation 
in the orthogonal direction can be neglected, 
with the consequence that the pressure is cons- 
tant in each cross section of the flow. This 
simplification has fundamental implications at 
physical and mathematical levels. In fact, the 
steady-state Navier-Stokes equations which form 
a coupled set of elliptic equations reduce 
to parabolic. The parabolicity of the boun- 
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SYNOPSIS 


Two dimensional flows where the boundary-layer approxi- 
mations are not acceptable are introduced and their specific 
mathematical and physical problems discussed, 

Inherent difficulties to the mumerical treatment of the 
partial differential-equations are exemplified and it is 
shcwn how increasing the Reynolds number imposes contra- 
dietory requirements regarding stability and accuracy, 

Purbulence models are introduced together with a numes 
rical solution procedure appropriate to these type of emos, as 
developed by |». B. Spalding et al. 
dary-layer equations has the physical effect 
that downstream values of any quantity do 
not affect the upstream ones. From compu- 
tation point of view this means that downstream 
values are obtained in a marching procedure from 
the upstream ones, starting from their given values 
at the initial station (initial condition) and on 
the boundaries (boundary conditions). (*) 

Though acceptable in many circumstances, as 
experimental values show, the simplifications are 
not possible in many important situations where 
the full set of two-dimensional elliptic equations 
must be used if meaningful results are sought. 
Some typical examples of separated flows or 
recirculating flows are show on fig. 1. In many 
cases, it is possible to divide the flow field in 
regions of full two-dimensional flow and bounda- 


(*) i.e. for velocity, usually a non-slip condition at the 
wall boundary and a finite value on the free stream are 
specified. 
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ry-layer type flow, as for example the impinging 
jet (Fig. 2) or the region of separated flow in a 
wall boundary-layer. (Fig. 1,d). 
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Such expressions as boundary-layer, sepa- 
rated flow, recirculating flow, . . . are sometimes 
confusing regarding their precise meaning. Thus 
we will classify two dimensional flows, after D.B. 
Spalding, in parabolic and elliptic flows according 
to the type of equations able to predict them, 
keeping in mind that, in this respect, parabolic 
flows are a simplification of elliptic flows. 

In this III part, elliptic flows will be considered 
only to stress present difficulties regarding their 
prediction, and to emphasize, by contrast, what 
simplifications were involved in parabolic flows. 
Detailed discussion regarding numerical solution 
procedures will not be attempted because it would 
take us too far out of the scope of the present 
course. 


1— THE MAIN EQUATIONS 


The Navier-Stokes equations, for an incom- 
pressible fluid without body forces when the flow is 
two-dimensional and steady (1-15) * reduce to 


dV od NV, 1 
V, | | a 1 
dx; UXa é OX 
BY 0" MV 
Rd e SS SR (1) 
Ux, 0X, 
d Vo O Vo 1 0 P* 
Dao ag a dm o 
| 0X, . 
02 V,  d2V, 
+uld + 
dX| dXo. 
« 


and the continuity equation to: 


d) V, d Vo 
+ — = 0 (2) 


dX| d Xy 


These equations form a coupled set of quasi- 
“linear elliptic equations, and we will call these 
flows elliptic flows (in the boundary-layer appro- 

d 2 
ximation ) ; ==0 in the first equation and the 
X 


* 


second is reduced to =0or P*=P*(x,); 


X9 


* [-15 means expression 15 of part 1; as a rule, Iorll 
preceeding a number refer to expressions in part I or II. 
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because of the simplifications, the equation 
become parabolic). 
If the flow is turbulent the equations are, for 


the mean values, in steady state, from (1-23) 


OU, UU, 10P o2U, oU, 
U— + U,—= — -» —S |—- 


9 2 4 
0X | 0X, 2 OX OX, UXg 
= ; Gs d 
-|- >u, + u, vs | 
dx, 2 0% 
oU, 9U, 
dx | 0X 
1 9P J? U, 0* U, 
ES — + y — - — 
o do dx? dx 
ds | 2 
O — — 
o E 4 
ae Regi 5) E dm u, (3) 
| 0X, UXs * / aa 
dU, oU, 
— +-=0 (4) 
0X, 0X, 


2— THE MATHEMATICAL PROBLEM 


2.1. — Boundary conditions 


As already refered, the fundamental equations 
in steady-state form a coupled set of quasi-linear 
elliptic equations, regardless of the flow beeing 
laminar or turbulent (*). The solution of this set 
implies the specification of the velocity or its 
derivatives on the whole boundary (**). Asan 
initial value problem (Cauchy problem) the pro- 
blem is not well posed as shown by Hadamard 
for the linear case. 

The need to specify the unknown on the whole 
boundary has the physical meaning that the value 
in each point influences the values in all others 
(which is typical of elliptic problems). From a 
prediction point of view, it means that in many 
important situations a boundary value is to be 
specified at a location where we would like to 
use the equations to get it by calculation. 


(*) In turbulent flow, we must stress, the instantaneous 
(time depend) equations are parabolic. They become 
elliptic only after the averaging process, and in the 
assumption of a steady mean flow — see appendix., 

(**) Mathematically, specifying only the derivatives 
would only insure uniqueness within a constant. 


TECNICA N.º 407 


As a simple example, which illustrates the 
practical difficulty, let us consider the fowin a tube. 
The elliptic nature of the problem imposes, for 
its solution, that an inlet and an outlet profile 
be specified, besides a null value on the tube 
wall. However, to specify the outlet profile is, 
in many cases, to specify what we would like 
to calculate. 

From a practical point of view, we will try to 
avoid the difficulty, specifying the outlet profile 
far downstream from the cross section where the 
velocity is to be predicted, on the assumption 
taht, ifthe boundary is far enough downstream the 
values specified there would not affect the predic- 
tion. This is not true from the very nature of 
elliptic equations, so an error is always involved. 


Fig. 2 


Another example is the impinging jet (Fig. 2). 
The mathematical problem implies the specifica- 
tion of the unknown on the boundary. Though 
the pratical way would be to specify in abc, def 
the values from a potencial flow solution, in cd, 
fg and ha the values from a boundary layer 
solution, we easily realize that such a procedure 
can entail significant errors (*). Another typically 
simple situation is depicted in Fig. 3. Values of 
velocity are to be specified in AA and BB. Of 
course, if the square cavity repeats itself 
periodically, the assumption of an equal profile 
in AA and BB can be adopted and its value 
found by iteration, starting from a tentative 
profile. 

The practical difficulty just mentioned for simple 
cases (where some remedies have also been 
devised) shows that even if the solution of the basic 


(*) We omit reference to part of the boundary which 
coincides with a solid wall. There, the non slip condition 
is obvious. 
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equations were known, a better understanding 
regarding the way of specifying boundary condi- 
tions would be still necessary. 


2.2 Numerical schemes 


Assuming that we know how to specify the 
boundary values, we must still face the problem 
of solving mumerically the set of equations. In 
this respect, a long way is to be followed before 
reliable and cheap procedures are found (despite 
many claims that these goals have already been 
achieved). 

Regarding numerical procedures, known me- 
thods have not yet been able to overcome the 
difficulty of a predictable effect on accuracy- 
-stability-computer time of increasing Reynolds 
number. To clarify the nature of this problem 
let us consider a simple example: the approxi- 
mation of convective terms by finite differences. 
The example is based on a model equation and 
follows closely the discussion of Fromm [3]. 


Consider the equation 


Í ;2W 2w 
es rue mv( 5 + o ; ) (5) 
dt 0X, dx] dx, 


where w is any convected property (for example 
the vorticity), » is a constant (diffusivity) and Uis 
the fluid velocity which will be assumed constant 
for the sake of the example. 

From a Taylor's series expansion around the 
point x =idAx, t=nÃt, we get: 


n 
n n dw 
0X, / 1 


Piada (e) 
5 ae 
2 0x) 
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A +? , n (6) 
Mi À. à PA 
2 Not” (ij 


(where the subscript i represents the spatial 
coordinate and the superscript n the time where 
the derivatives are to be calculated) 

And the left hand side of (5) is approxima- 
ted by (*) 


n-+I n 
dw dw o Ecos não 


D+ 47) 


ot dx At 


2 
sy MW be d w 


AX, 2 àx, / dO x7 

If we take the lowest order approximation, i. 
e. only the first two terms in the r.h.s. of (7) 
the error 15 


A x, / A à a” 
= [= nd e 


which shows that the approximation introduces 
an artificial diffusion of the above value. 

If we now consider the right hand side of (5) 
we see, regarding the errors involved in its 
finite-difference approximation, that our finite- 
difference scheme is the analogue of a different 
equation which is 


dW OW 
dt 0x1 (8) 
UAx A : 
= (+ am Sede 
2 N à x, 0x4 


though, the limit 4x, At — O shows the finite- 
-difference scheme to be consistent. Practically, 


A x, Udt 
Ax, Atare finiteandunlessU— |[1——— I<<y 
1, 2 A X| 


we shall expect important errors in the calculation. 


(*) on the assumption of v= o 
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Of course, we can choose 1 — a = 0, but 
this is possible only in our simple example 
where U is constant. 

In the general case U is not constant and can 
assume both positive and negative values. If w is 
the vorticity, and the equations are written in 
non-dimensional form, the non-dimensional diffu- 
sivity is 1/Re, so if Re increases, At and 4x, must 
decrease if accuracy is to be preserved. On the 
other hand, regarding stability, the total diffusion 
coefficient 


as At 
pod A (1— U =) 


2 8x, 
must be positive. If not, the basic solution would 
diverge. 

If this coefficient is approximately zero, we 
have the typical problem of hyperbolic equations 
which are usually stabilized by introducing an 
artificial viscosity. 

Summing up: taking into account only the 
finite-difference approximation of the convective 
terms, an error is involved due to the finite-size 
of the mesh, i.e. 4x, 4t. This error produces the 
equivalent of a diffusion, with a diffusivity 
coefficient of value 


U dxu (1 


2 à xa 


U àt ) 

In order to get accuracy, this coefficient must 
be small compared to the true physical diffu- 
sivity. This means that, with increasing Reynolds 
number the numerical diffusivity must be reduced 
to preserve accuracy. However, numerical stability 
imposes a finite positive viscosity, regardless of 
it being true physical viscosity or numerical. 
Because true physical viscosity becames vani- 
shingly small with increasing Re, the only way 
out to stability is to accept numerical viscosity, 
at the price of a loss in accuracy. 

Though very crude (**), the example shows 
the conflicting demands on the finite difference 
approximations of this type of equations: for 


(**) Ata first sigth the reader would consider this an 
artificial example of elliptic flows because the model 
equation is time dependent. In fact, because of non-linea- 
rity, the solution of the final algebraic system may 
be found by iteration. This iteration is the equivalent of a 
time dependent term, which justifies the model considered. 
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stability, a finite, positive, viscosity must exist. 
If the physical viscosity » is not large enough, 
an artificial viscosity is needed which in turn 
strongly affects the accuracy, 


1 
Because» < > =" the strong influence of the 
e 


Reynolds number in the computational problem 
is easily understood. 

The above remarks apply equally well to lami- 
nar and turbulent flows, provided the turbulent 
viscosity or its equivalent is taken into account. 


3. THE PHYSICAL PROBLEM IN TURBULENT 
FLOWS 


The main mathematical difficulties have been 
briefly outlined, on the assumption that the Rey- 
nolds stresses are known. The law which these 
stresses must follow is the physical problem, 
which is far from being solved. 

By reasons often refered, the computation of 
the Reynolds stresses from the N.5S. equations 
seems presently an impossible task. In purely 
time dependent laminar flows the numerical sta- 
bility is lost, well before the critical Reynolds 
number is attained, unless an artificial viscosity 
is introduced which renders meaningless the com- 
puted values. In turbulent flows, the mean flow 
equations are obtained by the procedure already 
discussed in part I, introducing the additional unk- 
nown of the Reynolds stresses. These are to be 
postulated in a way similar to that discussed in 
part II. In this case, however, the turbulence 
model is necessarily more complex. Taking, for 
instance, the production terms in the turbulent 
kinetic energy equation we have now four, ins- 
tead of the two terms which appeared in the 
boundary-layer case. 

Besides, in the boundary-layer case there are 
a considerable amount of experimental values 
to guide in the formulation of models or to test 
them when, in the elliptic case, data is very limi- 
ted and unreliable. 

Because the physical situation is more complex, 
more data is needed to account for the larger 
number of influencing parameters. These, in turn, 
render more difficult the empirical correlations and 
stresses the need for more refined theoretical 
models and computational procedures which 
would ease the need for experiment. 
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Before we present the example of a remarkable 
attempt in such direction, we mention two diffi- 
culties : present techniques of hot wire anemome- 
try rely heavily, for its interpretation, on previous 
knowledge of the mean flow direction (and 
apply only in low turbulence levels). This 
means that (though reliable in boundary-layer 
flows) the basic theory of interpretation of hot 
wire signals must be revised to deal with 
elliptic flows (this is the subject of following 
lectures). Use of turbulence models, like a 
mixing length, which has been sucessful in 
many practical situations where a typical length 
exists (such as the distance from the wall) is 
hurt with the difficulty that in general elliptic 
flows no such typical geometric length is obvious, 
which renders adventurous to extrapolate results 
from boundary-layer flows. 

As a conclusion, we can say that a tremendous 
need exists for experimental data which would 
inspire and test turbulence models, and that 
improved turbulence models are needed to guide 
experimentation. 


4. A PREDICTION PROCEDURE 


4.1. Having in mind the difficulties refered 
previously we would like to briefly present the 
approach of Spalding et al. [4] develloped at 
Imperial College to deal with the prediction of 
elliptic flows. 

We will refer mainly to the thesis work of 
Wolfsthein [5] because the theoretical predictions 
were partially compared with experimental values 
by Russel in Manchester. 

Assuming a Boussinesq form, the Reynolds 
stresses 
(1, = 1,2) 


Ryj=tu;u, 


9U, oU, 
==" EM ) (9) 


dx j d X; |) 
and introducing an effective viscosity 


Peg = + He 


“* Obviously if wu, =0, all the expressions apply to 
laminar flows. 
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The basic equations (3) can be written in the 
form: 


TOW jP 
O aU, — 
2 o X| já d X| 
dd dd U, 
o ) efi ) Ei 
ia Ux, 
d d U, 
e cd mel B— — w) (10) 
7, X3 J Xj 
1 df vº d Pp 
on: = 1 U, w + —= 
2 )Xo O X 
e) ( ) U,+ 
= 2 o ) 
', À efi J X, 
d ') 
EE “ef [a sEem) 
d | x. 
2 
where 
dU, J U, 
W = o (11) 
d X, J Xa 
is the vortici!y, and 
V=U,+U; (12) 


is the velocity squared. 

Taking the derivative of equation (10) in x, 
and subtracting the second after derivation in x, 
the pressure term can be eliminated, and equa- 
tions (10) are transformed into: 


9 
dw dwW Jr 
EU, Fey = — e W) + 5 
dxX| dX, dx, dx; 
(13) 
where 
SU J? pa 
) “ eff 
Se=2(2—+—— 


0X , 0X| 0X, 


9oU dj 2 uu JU )2 it 
Ea nt et e eti ) (14) 


E) 
OX Ux UX| dx ' 


Introducing now the stream function: 


FW = fou, d, — [e U, dx, (15) 
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The continuity equation (4) is transformed into 
2 

dx; dx 

J 


(16) 


Ss e 0 W 


Equations (13) and (16) are the two fundamental 
equations which determine the two unknowns: w 
and 'F, 

After w and W are known, U, and U, are 
found from (11) and (15) and the pressure P 
from (10). 

The basic equations (13), and (16) are of the elli- 
ptic type. As already refered, its solution imposes 
a specification of the unknown or its derivatives 
on the whole boundary. On any surface, S 


(LO =Evpi() 
.n/) 


de 


(17) 


is the boundary condition for 'F. 
For w we get, from the continuity equation 


2 W 
[E ) e eW (22) 
s 


If S coincides with an impervious wall (non 
slip and zero normal component of velocity) 
equation (17) reduces to 


(18) 


o F. asi 
) = 0 or 'É, = const. 
dn /w 

Remark that the boundary conditions are 
coupled what can give rise to some additional 
numerical difficulties). 

4.2. For the numerical solution, the equations 
are put into finite difference form which are 
derived from the integration around the given 
point and assuming a linear variation of the 
unknowns in the domain —the original reference 
[4] must be consulted for details. 

The non-linear, algebraic set of equations 
which has to be solved can also be obtained 
througth formal finite -difference techniques using 


CE qm EN 


fined within an arbitrary constant as is well known. 
Relation (17) implies, of course, that U, and Us are 
specified on 5. 

(**) This relation is just a consequence of (16) with x, 
taken along the normal to S. As such, is a general rela- 
tion which applies in any surface. This gives the general 
boundary condition for w. Besides, if the velocity profile 
is given along S, this defines 'W', from (15) within an arbi- 
trary constant. 'W is the boundary condition for !t. 


) from (15). W is de- 
/ 8 
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uncentered differences. The use of uncentred 
differences (which the authors justify on physical 
grounds: a point in the field is more strongly 
aífected by values upstream than downstream) 
allows the stabilization of the numerical scheme, 
though the formal proof of such range of 
convergence is not known. As should be expected, 
from our previous example, the scheme introduces 
an artificial viscosity which gives rise te numerical 
diffusion [6] (*). 

However, provided the limitations are clearly 
understood, the method is very interesting and 
usefull, 

4.3. In turbulent flows, the Reynolds stresses 
are found using the assumption that all the 
quantities which appear in the turbulent-kinetic- 
-energy-equation are expressible in terms of the 
turbulent-kinetic-energy (k = q”) and a length- 
-scale., 

By dimensional analysis [5], the following 
expressions are found : 


Viscous dissipation: 


2 


| 
Ju,' k * 
ti - A — e. 
fra Cp hs 


(19) 
Diffusion : 
à (1 a ) o | u, dk | 
pi dx, 2 Piu dr” Me dx; d, 0X; 

(20) 

Reynolds stresses: 

d u, d a) 

fai =— READ — 21 
eu; U, = 4 Es à Fx (21) 


with », found using the Prandtl-Kolmogorov 
assumption: 


ta 4 (22) 


ty = CG Pk 


with these assumptions, the turbulent-kinetic- 
-energy equation becomes 


2 dk | / vt dk 
2 (Uix — dx Lud ii 


(23) 


2 JU JU, JU k!? 
i Ra i 
2 BG + Ta) 7 + Cp 


(*) Which, we think, explains its stability. 


After Spalding, if the Reynolds number of 
TA 


4 


v 


turbulence defined by R, = is very high, 


the length scales lpr ly, are related to y (distance 
to the wall) by 


E ul 


— == — = const = 1 (24) 
y IL 
which are asymptotic forms of 
— =1— exp (— A,.R,) 
y 
Ip 
E (25) 


On the above forms, all the expressions which 
define the turbulence model depend only on 
unknown constants which are hopefully univer- 
sal and can be found by experiment. 

Such a first tentative set has been proposed 
by Wolfshtein as: 


A, = 0.016; Ap = 0.263; C, = 0.22; 
Co = 0.416; c, = 1.53 


The scarcity of experimental data does not 
allow a full comparaison between predicted and 
measured values, so a final judgement is impos- 
sible. However, for the impinging jet, Russel 
has tested the predictions of Wolfshtein and the 
results are remarkably good if proper account 
is taken, both of the complexities of the model 
and of experimental verification. As a rule, mean 
velocity profiles are fairly-well predicted ; turbu- 
lence-kinetic energy profiles show, sometimes, 
opposite trends. 

Consideration of geometries other than the 
impinging jet are crucial for the test of «univer- 
sality» of the constants once they have been 
found from the impinging jet data. 


4.4 The above model still have a strong 
practical limitation because the parameter y is 
sometimes of dubious specification. 

To avoid this difficulty, Spalding proposed an 
improved model where the length scale is got 
through an additional differential equation. 

Spalding introduces a new unknown, the mean 
turbulent vorticity squared: 
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and a new equation to be verified: 


à Q “1 / v, ' ad OQ 

se E] 

=| +5)|+ 

dx; dx, | A 0) dx, 
dw? Eu ai 

+Ciu (—) +c,wvV0-c,0% 
0X, 


This equation is proposed, not deduced. Toge- 
ther with the additional assumption 


aca PE. 
Va 


Cy EE a,s[l-exp(— A, R |) | 


results 


Cp=apll-exp(— AGR) ] 


k!2| 


t y 


with the values for the constants 
0, — PRE ;, gg = 0.227 : A q == 0.016 : n= 1 


a,=0.416;A,=0.263;C,=2;€C,=1;C,=1 


The constants are suposed universal. The above 
values must be considered tentative. 

The new model was tested in boundary-layer 
flows with promising results. 


5. CONCLUSION 


Compared with boundary-layer flows, elliptic 
flows embrace a much wider class of important 
practical problems. However, their complexities 
are much greater and their sucessful understanding 
and prediction have to overcome experimental, 
physical and mathematical difficulties. Some of 
them were refered. 

Tentative turbulence models and a mumerical 
scheme able to reasonably predict known expe- 
rimental data already exist, but extensive devel- 
opment is still needed. 
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APPENDIX 


Classification of Navier-Stokes equations and 
their simplified forms. 


Because contradictory statements appear in the 
literature regarding the classification of the N. S.- 
equations it seems useful to add a short note on the 
subject. 

A formal proof of the type (*) being out of the scope 
of this course, we restrict ourself to one of the momen- 
tum-equations in the incompressible case (which is 
of second order) as follows: 


1 oP O 


dt PEN 2 : º 
dX; dXj p OX; 


d P ua 
If, for the moment, we assume that PY whenelimi- 
O &1 


nated throught the continuity equation does not intro- 
duce second order derivatives in x,, equation (a) is of 
the parabolictype in t because all the spacial coordina- 
tes appear in second order derivatives with the same 
sign and the time t only appears in a first order 
derivative whose coeflicient has the opposite sign to 
12 
the one ado, 
a 
In this case, the solution is obtained marching in 
time i.e., «what happens in the future does not affect 
the past». 
Equation (a) is the model ofthe full, time-dependent, 
incompressible, N.S. equations which are parabolic. 


dV 


IH the motion is steady, '—o and (a) reduces to 


ot 


d V, Pk Vs 1 OP 
(b) E + — —— = 0 
dXj; 0x, 5 DX j 


In this case, all the independent variables appear 
in second order derivatives with coefficients of the 
same sign: the equation is elliptic. We remark that, 


dV 
if we introduce a time average in (a) such that '=o 


dt 


we get the same conclusion. 

If we now consider the boundary-laver approxima- 
tion of (a), the second order derivatives in x, are drop- 
ped and (a) reduces to 


dV, JV dV, 
Eraldo Rh: Semiliado RR as 
2 
dt 0X Xg 
(c) 
9º V, 1 9P 
0X, - 0X| 


(*) We recall that the classification of partial dif- 
fential equations is based on the higher order deriva- 
tives which apper in each independent variable. 
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In this approximation, only x, appears in second 
order derivatives, and t and x, only in the first order 
ones, 

Equation (c) is parabolic in the broad sense. 


dv, 
In steady state, = 0 and we have 
0Vi ay OV 2v P 
é dei Tá “o 7 e! 
l ( 2 d Xo p UX| 


and the equation is parabolic in x, in the restricted 
sense or simply parabolic (the physical meaning of 
the equation implies V, > o which is also needed if 
equation (d) is to have a bounded solution). 

Because (d) is parabolic in x,, which isa space 
coordinate, the implication is that downstream values 
do not affect upstream ones (cf. with the role of time 
in equation (a)). 

In the general, compressible case, elimination of P 
throught the continuitv equation can introduce higher 
order derivatives whose coefficients can change sign, 
and accordingly the equations can change type. A 
formal and complete discussion of this topic is not 
known to us in the general case of viscous flows. In 
inviscid flows, the equations change to hyperbolic at 
super-sonic velocities, being of mixed type on the 
sonic-transition. 


NOMENCLATURE 


A,C — constants 


k ——turbulent kinetic energy 

| —length scale 

P | — pressure 

U, —mean velocity component in direction i 
u, — fluctuating velocity component in direction i 


V, — velocity component in direction i 
w  — vorticity 

x;, x; — orthogonal coordinates 

y —— distance to the wall 


u — — dynamic viscosity 

“, —turbulent dynamic viscosity 

y — — Kinematic viscosity 

v, —turbulent kinematic viscosity 

2 — specific mass 

3 + — constant 

S —boundary 

W — stream function 

O —turbulent vorticity squared (mean) 


Obs: — an overbar means always a time average 


355 


REFERENCES 


1] —J. J. D. DOMINGOS—The fundamental equations 
of turbulent incompressible flows, Tecnica 406, 
p. 287, 1971. 

(2! — J. H. WHITELAW — Prediction procedures for 
two-dimensional boundary-layer flows: equations 
and assumptions, Técnica 406, p. 299, 1971. 

3] — J. E. FROMM — Practical Investigation of Convec- 
tive Difference Approximations of Reduced Dis- 
persion. 

The Physics of Fluids Supplement II, 1-3, 1969. 


[4] — A. D. GOSMAN, W. M.PUN, A. K. RUNCHAL, 


D B. SPALDING and M. WOLFSHTEIN — Heat 
and Mass Transfer in Recirculating Flows — Aca- 
demic Press, 1969. 


[5 — M. WOLFSHTEIN — Convection processes in tur- 


bulent impinging— jets Imperial College of Science 
and Technology — Department of Mechanical En- 
gineering — Report SF/R/a2. 


[6] — A. K. RUNCHAL, D. B. SPALDING, M. WOL.- 


FSHTEIN — Numerical solution of the elliptic 
equations for transport of vorticity, heat and matter 
in two-dimensional flow — The Physics of Fluids 
— Supplement II, 1 21-1969. 


FEIRAS, CONFERÊNCIAS, EXPOSIÇÕES E CURSOS, NO MUNDO 


OBJECTIVO E PROGRAMA 


8emes Journées Européennes de la Commission Inter- 
nationale des Irrigations et du Drainage 


Recherches Européennes Sur La Résistance a L'incendie 


Técnicas de Entrevista 


= ——— 


4eme (Conference Internationale sur la Physique des 
Plasmas et la Recherche Concernant la Fusion Nucléaire 
Controlée 


Achema 71 — 17.º Congresso e Exposição de Aparelha- 
gem Química 


Small Gas Turbines for Helicopters and Surface Transport 


IMEBORON — International Meeting in Boron Com- 
pounds 


vil Assembleia Geral do CIB e V Congresso do CIB 


356 


LOCAL, DATAS, ORGANIZAÇÃO, 
INFORMAÇÕES 


Aix-En-Provence, 14-19 Junho-71 
Secrétariat de L'A.F.E.I.D. 19, Av. du 
Maine, 75-Paris, 15, França 


Paris, 15 Junho-71 
ITBTP 9, Rue la Perouse, Paris 16 


Lisboa, 16-18 Junho-71 
INI 


Madison, Wis., 17-23 Junho-71 


Frankfurt, 17-24 Junho-71 
Câmara de Comércio e Indústria Luso- 
Alemã 


Paris, 21-22 Junho-71 

Londres, 24-25 Junho-71 

Montreal, 28-29 Junho-71 

Pasadena, 1-2 Julho-71 

Direcção do Serviço de Materiais da Força 
Aérea Portuguesa 


CR Eme E “cs 


Castle Libice, 21-25 Junho-71 

Institute of Inorganic Synthesis, Czechos- 
lovak Academy of Sciences, Rez Near, 
Prague, Czechoslovak 


Paris, 21-30 Junho-71 
CIB — Postbus 299-Weena 700, Rotter- 
dam, Holanda 


TECNICA Nº 47 
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IV— AN INTRODUCTION TO HOT WIRE ANEMOMETRY 


RESUMO 


Justifica-se o uso de fios quentes em anemometria, em 
particular descreve-se o método de calibração destes, e 
interpreta-se a resposta do anemómetro de fio quente em 
escoamentos de baixa turbulência, 


1. INTRODUCTION 


1.1 — In fluid dynamics, measurements are 
needed in order: 


a-—to satisfy the requirements for further 
data resulting from the study of compli- 
cated flow situations. Generally, one starts 
with measurements assuming simple models 
for the flow, and then these models may 
be developed in such a way as to represent 
better the situation under study ; 


b—to provide confirmation of theoretical 
models. In this case, the scientist has to 
reproduce the assumptions of the theory, 
and the measurements will become much 


harder. 


1.2 — In the flow of fluids, the velocity at each 
point may be assumed to be the sum of two 
components, one being the time mean velocity 
and the other being the fluctuating velocity 
whose time mean value is zero. Choosing a 
coordinate system (O,x,,x,,xs) as indicated in 
figure 1. The instantaneous velocity at point O 
is then considered as the vectorial sum of the 
time mean velocity (U,,U.,U,) and the fluctua- 
ting velocity (u,,u,,U,). 


1.3— If the fluid flow is turbulent, the compo- 
nents of the vector (u,,u,,u.) are randomly 
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SYNOPSIS 


This paper provides an introduction to the physics of 
hot wire anemometry, The use of aires as tranaducers for 
turbulence measurements is justified, and for flows with low 
turbulence the means of signal interpretation is discussed, 

Techniques of calibration are aiso described, In parti- 
cular, the constant temperature (constant resistance) and 
constant current hot wire anemometers are presented, and 
their uses are commented on, 


distributed. When using statistical theories, 
turbulent flow is treated as an ergodic phenome- 


X3 


VELOCITY VECTOR AT TIME É, 
THROUGH POINT O 


(UU Us) 


Figure 1 - Components of velocity relative to the 
coordinate system (0,X,, Xy, Xsz) 


non, i. e. one assumes that time averages and 
ensemble averages are equal. 


1.4 — The study of turbulence as random 
phenomenon requires the use of transducers 
sensitive to the fluctuating component of velocity. 
The hot wire anemometer is one of the trans- 
ducers we may use to detect the fluctuating 
velocity in the flow. Other techniques do, however 
exist for making turbulence measurements, and 
in the lessons of this course, F. Durst provided 
also some information about laser Doppler- 
-shift anemometry applied to turbulence measu- 
rements. 
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1.5 — The information obtained from hot wire 
anemometer measurement can be processed to 
provide : 

a. — correlation measurements, in particular 
the components of the tensor u, u; intro- 
duced in the lessons of Prof. D. Domingos; 

b. — spectrum analysis and power spectra] 
density of the kinetic energy ; 

c. — r.m.s. turbulence intensity (*); and 

d. —the degree of isotropy and homogeneity 
of the turbulence field. 


2. THE USE OF WIRES AS TRANSDUCERS 


2.1 — The application of wires as transducers 
depends on the detection of the fluctuations in 
electrical resistance of the wires. The electrica 
resistance of a conductor depends on its dimen- 
sions and temperature, and the resistivity of the 
material. Table 1 lists those factors and names 
the coefficients used to specify the influence of 
the corresponding factor on the resistance. 
Table 2 shows the applications of resistance 
transducers corresponding to the factors on 
which the resistance depends. 

2.2 — The use of wires as anemometer trans- 
ducers is based on the relation between convection 
coefficients and the velocity of the fluid flowing 
around the wire. In hot wire anemometry, we 
use wires under conditions such that the resistance 
of the wire depends mainly on its temperature 
and thus on the heat transferred from or to the 
fluid by convection. 


The material of the wires used in anemometry 
should be such that: 


a. —it has a high temperature coefficient of 
resistance ; 


(*) — the rms. (root mean square) value of a function 
V (t) is defined as: 
12 


Pra JR A TP, 
9 = (V? (é) | | lim | vê (ya ) (1) 
' ' T “0 / 


—e 


In practice, for the purposes of measurement, the r.m.s 
value is calculated as: 


PN 1 1 “T, 1/2 
“meas | Vicio o)? = (+ | vê (t a) 48) 


N 


vêIdo 


where T, may be chosen sufficiently large that o is 


practically constant in time, and then: 


meas 


o 


“mes (3) 


358 


b. —it does not become oxidised easily when 
heated to high temperatures ; and 

c. — it is possible to manufacture small diameter 
wires whose characteristics are uniform along 
their length. 


2.3 — Table 3 records values of the temperature 
coefficient of resistance for some materials. Of 
these materials, the most commonly used in 
anemometry are tungsten and platinum. In spite 
of the temperature coefficient of resistance of 
copper being of the order of those of tungsten 
and platinum, copper is not used because it 
becomes oxidised even at low temperatures. 


3. HEAT TRANSFER FROM A HOT WIRE 


3.1—Let us consider an infinite wire, heated 
by the Joule effect, and let us suppose that heat 
transfer from the wire to the fluid exists by 
convection only. Due to the cylindrical sym- 
metry of the system, it is obvious that the 
temperature of the wire should only depend on 
the component of the flow normal to the wire. 


3.2 — An energy balance applied to a control 
volume containing the hot wire gives: 


Q,+0,=0,+0,+Q,+Q,, (4) 
where : 


Q | — heat generated by Joule effect 

Q, = heat generated by Thomsom effect (due 
to the existence of temperature gradients in the 
wire). This term is of the type: 


dT 
Qp= 8 die “4 , (5) 
dx 


where 7, is a constant. However this term is so 
small compared with Q, that it may be neglected, 

Q, = heat absorbed by the wire to change its 
temperature. It depends on the temperature 
variation of the fluid as a function of time, and 
is important only when using the hot-wire in 
non-stationary conditions. 

Q. = heat conducted along the wire to the 
hot wire supports. Hot wire probes cannot be 
manufactured with infinite wires since they have 
to be supported. Also we need to make local 
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méasurements, and so the wires must be short 
(up to about 2 mm). Thus instead of constant 
temperature along the wire, the temperatures at 
the extremities of the wire are the same as the 
temperatures of the supports (see figures 2 and 3). 


SUPPORTS OF WIRE 


Ds OR GOLD PIATED wIRE 


Figure 2— The probe for hot wire anemometry 


Q, = heat transferred to the fluid by convec- 
tion. The heat transfer ccefficient depends on the 
velocity of the fluid, i. e. on Reynolds number, 
and for values of Re < 0.5 natural convection 
predominates over forced convection. If we want 
to measure fluctuating components of velocity in 
turbulent flow need to use the smallest possible 
diameter for the wire, in order to obtain a high 
frequency response; with tungsten wires this 
diameter is about 5 m, Which gives us a lower 
velocity limit of about 1 m/s. 


HOT WIRE TEMPERATURE 


Q, = heat radiated from the hot wire to the 
outside medium. This may be written: 


Q =50, (e T* — à To!) , (6) 


where : 


S = external surface area of the hot wire, 

o, == Stefan-Boltzmann constant, 

e = superficial emissivity of the wire at the 
temperature T (Ok), and 

« = absorptivity of the surrounding medium 
at temperature T, (Ok). 

This term may also be neglected in comparison 
with Q E 

Thus equation (4) may be written: 


Q,=0, + Q, + Q (7) 
Considering the relations that define Q y e Bee M 
and Q, per unit length of the wire, we have: 


Jº 


PT. 
A 


od T 1) ( o T 
q T “ay by 


+rdh(T-To) (8) 
If we assume that the process is stationary, 
then Q, is also negligble. 
Taking in account the boundary values, equa- 
tion (8) can be solved in order to obtain the 
temperature along the wire in function of time. 


NORMAL FlOw TO HOT WIRE 


— - —- YAWED FLOW TO HOT WIRE 


Figure 3 — Distribution of temperatures along the hot wire 
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3.3 — The probes used in hot wire anemometry 
consist of two supports across which may be 
soldered the wire whose dimensions are generally: 


d—from 3 to 74, more frequently 5 £; 

|— about 2mm, but in general 1 > 300 d. 

In figure 2 is sketched the characteristic form 
of a hot wire. 

To define the length 1 of the wire, to reduce 
the influence of the progs on the velocity 
measurements, and in some cases to facilitate 
the soldering of the wire to the supports, is gold 
or copper-plated at the ends of the sensitive 
portion. 


4. HOT WIRE CALIBRATION 


4.1 —Hot wires are calibrated in stationary 
conditions for which equation (7) may be written: 


Q,=0,+ OQ. (9) 


Neglecting the heat transferred along the wire 
to the supports (Q.), we may write: 


QO,=Qh, (10) 
OT: 
RIÉ=kh(T—T,). (11) 


A linear relationship can be assumed for the 
temperature dependence of the resistance R (T): 


R=R,[1+T(T—T), (12) 


where 'V is the temperature coefficient of resis- 
tivity of the material of the wire. 


Calibration curves can be obtained by placing 
the hot wire in a flow with constant velocity U, 
which gives us a constant value for h. Equations 
(11) and (12) show that there are then two 
variables R and I, we may also consider the vol- 
tage E between the extremities of the hot wire: 


E=RI. (13) 
Then we have calibration curves in spaces (R,1), 
(E,R), and (1,E) at constant velocity whose shapes, 


shown in figures 4,5, and 6 are given by equations 
(11) and (12). 
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u, 
R(n) U 
RESISTANCE MAXIMUM 
(Tempero ure ã Po INCREASING U ( A ) 
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! AT 
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PP Qu U= 
R, = 
ns E 
[e] 


Consrant 


WORKING AT I CONSTANT (with saturation j 


(wi th satwuratiom) 


IT (ma) 
R(1) = KA R 


KA- x Re 
Figure 4 — Calibration curve R (1) 


dá INCREASING U,À 


i UY WORKING AT 
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Vo 
A WORKING AT CONSTANT 
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72 necã SPENNES Ar CONSTANT 
; vaTaAGE E (with saluralion) 
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é E (Volt) 
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e dd q 


Figure 5 — Calibration curve I (E) 


I(E) = 


À 
psi Vs=o 4 y, 
INCREASING U, À 
E as 
apa is 
/ E ge tdo 
U= 00 


WORKING AT 
E constant (with salurct; am) 


E (volt) 
Figure 6 — Calibration curve R (E) 


4.2 — When we are using a hot wire for 
turbulence measurements, these calibration curves 
are not needed as resistance and current will 
not be fluctuating independently during the 
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measurements. Hot wire anemometers are built 
in such a way as to keep one of the variables R 
or I constant during operation; for such units, 
the calibration curve will simply be the other 
variable (I or R respectively) as a function of 
the velocity, or the voltage E as a function of 
the velocity E (U). 


4,3 — In figures 4,5 and 6 are shown sections 
across the curves corresponding to: 

a. — hot wires working at constant resistance which, 
as can easily be seen from equation 23), is the 
same as constant temperature ; 

b. —hot wires working at constant current; 

c.—hot wires working at constant voltage. 

We can conceive an infinity of sections corres- 
ponding to hot wires operating under many 
other conditions; however, the problem is to 
guarantee that we are working in these condi- 
tions. 

The most developed methods are the constant 
resistance (constant temperature) and constant 
current ones; these will be discussed in the next 
paragraphs. 


OUTPUT 
E (volt) 


Figure 7 — Constant current hot wire anemometer 


5. CONSTANT CURRENT HOT WIRE ANE- 
MOMETER 


5.1 — To obtain a hot wire working at constant, 
current, the wire is generally mounted in series 
with a large resistance so that the resistance 
fluctuations of the hot wire may be neglected 
compared with the series resistance. Figure 7 
shows the circuit of the constant current hot 
wire anemometer. Using the nomenclature in the 
figure, the current through the hot wire is 


E, 


| 


(14) 
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and a fluctuation of R gives a fluctuation of 
current 


E 


ET cu seg nei 
(R + Rs + Rg? 


AR. (15) 


If Rç is so much bigger than R, that AR« 
«R + Rg + Rg, then ÁIDO, and we say that the 
hot wire is working with constant current. 


5.2 — From figures 4,5, or 6 it is seen that hot 
wires working at constant current are saturated 
at high velocities, i.e. at high velocity the 
difference in the output from the hot wire unit 
for a given change in velocity is quite small, so 
that reproducibility in the measurements is hard 
to obtain. Also at high velocities, since the hot 
wire temperature is lower, the thermal capacitance 
of the wires becomes more important. This 
problem can be partly overcome by working 
with high currents at these high velocities; but 
a problem then arises that if there are large 
velocity fluctuations, or if the wire is removed 
from the flow without reducing the current, the 
wire may easily be oxidised (and broken). 


6. CONSTANT TEMPERATURE HOT WIRE 
ANEMOMETER 


6.1 — To obtain a hot wire working at constant 
temperature, it is necessary to use a feedback 
control system to supply a variable voltage to 
the wire, in such a may that its resistance is 
held constant. This can be acomplished with the 
circuit shown in figure 8 which operates on the 
following basis. 


OuTPUT 


| E (volt) 


As — AMPLIFIER WITH GAIN G 


Figure 8 — Constant temperature hot wire anemometer 


The hot wire is connected as one arm of a 
Wheatstone bridge, and a fluctuation in resistance 
of the wire causes a fluctuation in the voltage, e, 
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between A and C. This voltage when amplified 
by amplifier Aq is used to change the voltage 
between B and D in order to obtain e= O again. 

Assume that when the bridge is balanced, i. e. 
when e=o and R =R,, the output from the 
amplifier is E,. The voltage fluctuation e due to 
a fluctuation of velocity may be written. 


fi=R. 


e=--E, (16) 
2(R+R,) 


and if G is the gain of the amplifier, the output 


voltage is 
E=E,+G.e. (17) 


Combining (16) and (17), the hot wire resis- 
tance is given by 


E, 
cra(1- 0) 


Es + 18) 
R, € (1 E, 
— 2 = í 


With a very high gain amplifier, we have 


lim e” ; 
R, (19) 
G > oo 


and the hot wire is working at constant resistance, 
or constant temperature. In practice, the gain is 
limited by the noise from the circuit which may 
introduce alterations in the output from the hot 
wire unit. 


6.2 — At constant resistance, the calibration 
curves may conveniently be expressed in the form 


E =f(U) , (20) 


assuming flow normal to the hot wire. 

From dimensional analysis for the heat transfer 
from a hot wire, assuming only forced convection, 
we may write: 


N,=F(P,R,ç) . (21) 
With : 


N, = Nusselt number 
P, = Prandtl number 


R . = Reynolds number 


defined for the fluid properties and the diameter 
of the wire. 
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King's law for the relation (21) is of the type: 
N,=A+BVR, , (22) 


where A and B depend on Prandtl number, 
but may be assumed to be constants in the case 
of air (for which the Prandtl number may be 
taken as uniform). A more accurate represen- 
tation of the results compared with experimental 
measurements is given by Collis and Williams's 
law 


OP fi 
( 2t, 


— 17 os 
) =A+B,R, , (23) 


where T and T, are respectively the temperature 
of the hot wire and the temperature of the flow. 

As the power converted to heat in the wire 
working at constant resistance is proporcional to 
E,, we may write: 


E=ES+HBU" , (24) 


for flow mormal to the wire. This law comes 
from an overall heat balance and from King's 
law, and the coefficients B and n depend on the 
method of calibration. E, is the output voltage 
at zero velocity. 

King assigned to n the constant value 0.5, 
while for air Collis and Williams used 0.45 and 
stated that this value depends on velocity. 
Calibrations made by Davies and Bruun, and 
Kjellstrôm and Durst showed that n depends on 
the velocity U of the flow, decreasing as the 
velocity increases. In their calculations, n was 
found using the outputs E, and E, from the hot 
wire at two distinct velocities U, and U,; then, 
since 


E, =E BU) (25) 
and 
E, = ES + BU; (26) 
we obtain 
log (E? — E2) — log (E2 - E?) 
aa. É ai: > 


log U, — log U, 


6.3 — In the case of wires yawed relative 
to the flow, the response of infinite wires depends 
only on the component of velocity normal to the 
wire. Thus equation (24) will be written: 


EL =E2+B(U.cosa)”, (28) 
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where «a is the angle defined by the direction 
of the flow relative to the normal to the hot wire, 
and E, is the output from the hot vire anemo- 
meter when the probe is placed in a flow of 
velocity U at a yaw angle a. 

For wires of finite length, Champagne, Kjel- 
Istrôm, and others used the relationship 


E2 = E + BU" (cosa + Kºsin”a)"? (29) 


which reduces to equation (28) for K =. 

K depends on the value of n, and quite dif-. 
ferent value for K are obtained when using 
the n values of King, Collis and Williams, 
or Kjellstrôm and Hedberg. It has been shown 
that K varies linearly with the diameter/length 
ratio of the wire. 

Davies and Bruun use the equation 

Ei = E? + BU" cos” « (30) 
where n is the exponent defined with the flow 
normal to the wire and m is a value which 
depends on the diameter-length ratio of the wire. 
In order to compare the results from equations 
(29) and (30), these may be expressed in non- 
-dimensional form by comparison with equation 
(24) for flow mormal to the hot wire. 

The non-dimensional quantities : 


Ex y (E) +[-(22)] 


( cos? a + k?sin?a)o 


| (31) 
and 
E SÊ A E 
EV (Ef 2) [eos 


(32) 


can be compared when plotted as a function 
of a, 


7. SENSITIVITY COEFFICIENTS OF HOT 
WIRES 


7.1 — Let us supose that a hot wire is placed 


in a turbulent flow, mormal to the direction of 
the mean velocity vector (see figure 9). The ins- 
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& HOT WIRE 
Figure 9 
tantaneous effective velocity of the flow may be 
written 
ma: d 2 
Vr=VMUt+uP+u,tu?, (33) 
i.e.: 


É 2.A 2 
u + us + us 


us 
Ua=UV/a +28 + U? 


Considering the series development of the square 


root term: 
u 

Ua=U) ++ + | 
U 

(35) 


and neglecting all terms above first order we, 
may write 


(34) 


É ausitid ada 
utu+u; 


2U* 


Um U-u, (36) 


Thus the hot wire is sensitive only to the 
fluctuations of velocity in the direction of the 
mean velocity U, normal to the wire. 

Let us now assume that the mean velocity 
vector makes an angle ” with the normal to the 
hot wire (see figure 10). 


| HOT WIRE 


Figure 10º 


Then we may write: 
u,=º9U (37) 
u,=0"F (38) 
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As the output voltage from the hot wire is: 


then, for small fluctuations of E 
JE=e= o E JU + did 0W , (40) 
JU JH 
and from (37) and (38): 
e=a 46? i (41) 
U U 
where : 
E: 
a (42) 
and: 
: 09 E 
= (43) 


are respectively the coefficients of sensitivity 
parallel and normal to the direction of the mean 
velocity U. 


Using equation (29), and assuming n, k, and 
N 


n 
t (= a are independent of U 


2 2 Due à. Jwrt 

E = E. + BU" (cos PAKk sin Ch) , (44) 
we obtain : 
SE n(E'— E2) 


x = () nando “atá 
dU 2E 


(45) 


;.90E — uk I) (E —EgsinzY 


b TE 2o——. 46) 
04 2E (cos'V 4 kº sin !F) 
7.2 — Let us assume that we have a hot wire 
placed in a flow at an angle of yaw Y with an 
instantaneous output voltage. 


u 


2 
U , (47) 


4, 
e; e x, U + É, 


to which corresponds the r.m.s. value 


4 


dá Ma o Syd 
«t- 4 — > + 24, HM] 9 x (48) 
U U 


E é 
: u | 


U? 


If the wire is now turned until its normal 
makes an angle (— 4) with the mean velocity 
vector, the voltage output will be: 


u, U, 


"8 49) 
a ( 


364 


with a mean square value: 
uu, 
nd imãs (50) 
U 


este, (51) 


From equations (47) to (49) we can obtain the 
fluctuating components u, and u,, and the values 


of the tensor u, u,. 


7.3 — More generally, we can assume that the 
voltage between the ends of the hot wire depends 
not only on the velocity U and the angle v, but 
also on the temperature T of the fluid (assuming 
constant temperature hot wire operation). 


Then 
E=E(U,y,T) (52) 


If O is a small fluctuation of T, then the 
instantaneous output voltage will be 


fi 


= 4 q “a & (53 
raiado RS dE +7 é ) 
where : 
d E 
ua p És. 54 
/ E (54) 


is the sensitivity coefficient of the hot wire to 
the temperature of the fluid. 


From paragraph 2 we can write 
E? 

= Tg T), (55) 
where T, is the temperature of the hot wire, 
and h* represents the forced convection heat 


transfer coefficient depending on U and “. 
For small fluctuations “ of the fluid tempera- 


ture 
2 E Ú Lo 
t=|—|l=1-(1— 
| E | = O) 
or: 
f| 1) 1 / 6 2 
E = Td É as ra 
RE | Loo dog E | 


(57) 


and we can neglect s if,and only if, <<To-T. 
For liquids it is much more difficult to obtain 
0 << T,- T'than for gases because the hot wire 
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temperature needs to be lower to prevent bubble 
formation in the liquid near the wire, and also 
because it is difficult to obtain a uniform tempe- 
rature in the system. Thus in cases where the 
effect of temperature is important, the measure- 
ments obtained need to be corrected, Alternati- 
vely, the hot wire anenometer may be tempera- 
ture-compensated by using a second system 
sensitive only to the temperature of the liquid to 
correct the output from the hot wire. 


8. TABLES 


TABLE 1 


Electrical resistance of materials 


Factors on which resistance 
depends 


dimensions of wire, 


which depend 


on deformation 


at constant 
volume 


by hydrostatic 
compression 


temperature 


TABLE 2 


Corresponding 


coefficient 


sensitivity 
factor 


(gauge factor) 


pressure 
coefficient 


of resistance 


temperature 
coefficient 
of resistivity 


Applications of resistance transducers 


Factors on which the resistance 
of the transducer depends 


deformation 


temperature, 
which is 
altected by 
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MES 


at constant 
volume 


by hydrostatic 
compression 


heat transfer 
by conduction 


heat transfer 
by convection 


total heat 
transfer 


application 


extensometry 
transducers 


pressure 
transducers 


thermal con- 
ducetivity 

and thermal 
diffusivity 

transducers 


vacuum pres- 
sure 
transducers 


anemometry 
transducers 


temperature 
transducers 


TABLE 3 


Temperature coefficient of resistance for some 
materials (*) 


v = temperature 
= resistivity at coefficient 
OC (1090 cm.) of resistance 
(n.em. ºC —1) 


Material | “0 


Copper 1.) 4. X 10-35 
E cp == 4.9 2. >< 10-68 
Manganin 44 5. XX 10— 6 
Platinum E Io. 39 X 10-35 
pas 5.6 5. XX I0-5 


Note — resistivity at temperture T (ºC) assu- 
ming a linear dependence for q (T) is: 


(D=e [+47]. 


NOMENCLATURE 


O — point of measurement 

Ox, (i=1,2,3) — coordinate axes through point O 
U, (i= 1,2,3) — time mean components of velocity 
u, (i= 1,2,3) — fluctuating components of velocity 
v (t) — theoretical value of any function of time 


Vmeas (t) — measured value of v (t) 


v (t) — time mean value of v (t) 

R — electrical resistance 

I — electrical current 

E — electrical voltage 

| — length of wire 

d — diameter of wire 

a — sensitivity of wire to the component u, 

É — sensitivity of wire to the component u, 

7 — sensitivity of wire to the temperature of the 


fluid. 


(*) from: HANDBOOK OF CHEMISTRY AND 
PHYSICS 43rd Edition (1961-62) Chemical Rubber Pubs. 
Co. 
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SOBRE O CRITÉRIO LIAPUNOV 
NA TEORIA DA ESTABILIDADE DE SISTEMAS 


RESUMO 


Definem-se as noções de estabilidade e efectua-se a sua 
interpretação geométrica. Exemplifica-se o método directo 
eo critério Liapunov num sistema mecânico, Comparam-se 
os diversos métodos de formação de funções Liapunov em 
sistemas não lineares, 


1. INTRODUÇÃO 


No final do século passado Liapunov [1] 
estabeleceu dois métodos de estudo da estabilidade 
de sistemas. No «Primeiro Método» procurou 
definir as condições de estabilidade através da 
resolução da equação diferencial de estado do 
sistema. Mas dada a dificuldade nessa resolução 
Liapunov definiu o «Segundo Método» em que 
se aprecia a estabilidade sem solucionar a equação 
diferencial de estado, baseando-se num critério 
energético do sistema físico. 

As ideias contidas no «Segundo Método», 
também designado «Método Directo» por se 
ocupar directamente do sistema sem atender à 
resolução da equação de estado, constituem os 
alicerces dos métodos modernos de estabilidade 
de sistemas, nomeadamente no controle automático 
ou na análise do comportamento transitório de 
redes eléctricas. 

De facto estabelecendo-se a energia de um 
sistema dinâmico em função das suas variáveis 
de estado, se esta função do tempo t tende para 
zero quando t cresce significa que a energia se 
reduz monotonamente e portanto o movimento 
do sistema tende para o repouso, quer dizer, a 
origem das coordenadas de estado é um ponto 
de equilíbrio. A apreciação da estabilidade resulta 
assim da «função energia», subordinada às 
hipóteses indicadas (expressas nas variáveis de 
estado e tender para zero quando t tende para 
infinito). Se se conseguir porém definir uma 
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SYNOPSIS 


We define the conceptions of stability and show their 
geometric interpretatione. The direct method and the Lia- 
punov criterion are applied on a mechanical system. The 
methods of formation of Liapunov functions in non-linear 
systems ure compared 


função que, embora não representando a energia 
do sistema, satisfaça às condições enunciadas, 
poder-se-á analogamente julgar sobre a estabili- 
dade. A essa função chama-se «função Liapunov». 

Estas considerações dão ao «Segundo Método 
de Liapunov» um interesse puramente histórico, 
pois o problema consta agora na formação de 
uma função Liapunov que prove a estabilidade. 
Diversos métodos têm sido propostos, não se 
tendo contudo encontrado um esquema geral e 
prático. O recente critério de Popov por outro 
lado deu à questão uma importância prática 
secundária. 

Pretende-se neste artigo concretizar as ideias 
essenciais da teoria da estabilidade, enunciando 
os conceitos básicos dos métodos apresentados 
na literatura, e estabelecer uma sistematização 
didáctica do problema da estabilidade. 


2. DEFINIBILIDADE DE UMA FUNÇÃO 


Antes de se analisar o conceito de estabilidade 
é importante introduzir algumas noções básicas 
para o seu tratamento matemático. 
Chama-se norma ||x|| à distância || x || = (x; + 
1 


1 
= x) +... + x)" ou seja ||X || = (x! x)”, onde 
x" representa o vector transposto de X. 

Uma função continua V(x) diz-se definida 
(positiva ou negativa) na vizinhança da origem O 
quando nessa vizinhança possui derivadas parciais 
contínuas e em todos os seus pontos apresenta 


367 


o mesmo sinal (respectivamente, positivo ou 
negativo) excepto na origem onde é nula. 

a) Função definida positiva (apresenta um 
mínimo na origem): 


V(x) > 0 para todos 
V(0)=o0. 


x50 


b) Função definida 
máximo na origem): 


negativa (apresenta um 


V(x) <0 para todos xF0 


V(0)=o0. 


Uma função V (x) é semi-definida (positiva ou 
negativa) na vizinhança da origem O quando 
nessa vizinhança possuir derivadas parciais 
contínuas, for nula na origem e nos restantes 
pontos da vizinhança apresentar o mesmo sinal 
(respectivamente, positivo ou negativo) ou Oo 
valor nulo. 


c) Função semi-definida positiva : 


V(xg)—o 
V (0) =0. 


para todos x = O 


d) Função semi-definida negativa : 


Ví(x)<o 
V (0) = 0. 


para todos x == 0 


Uma função V (x) diz-se indefinida quando na 
vizinhança da origem O possui sinais positivo e 
negativo, isto é, quando não for semi-definida. 

e) Função indefinida : 


V E U, 
(x) < 


Por exemplo a função V (x) = x, + xá no 
espaço bi-dimensional R? é definida positiva pois 
só se anula na origem, para x, =x, = 0, e nos 
restantes pontos é sempre positiva. A mesma 
função no espaço tri-dimensional Rº é semi-defi- 
nida positiva pois se verifica o seu anula- 
mento para x, ==x, =x, =0 e além disso para 


sa comem fu - 
todos os pontos x, =x, = O com x; O. Aná- 
2 


logamente V(X) = X| + x) a x no espaço Rº 
é definida positiva e V (x) =— (x] + x5 east) 
no espaço Rº é definida negativa. A função 
V(X) = (x, + mr + x; no espaço Rº é semi-de- 
finida positiva porque se anula para x,=x,=x,==0 
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e ainda para x, == x, com x,;=0. A função 
V (x) =x, + x, é indefinida pois pode assumir 
valores positivos ou negativos conforme os sinais 
de x, e x,. 

Diz-se que a função V (x) assume a forma 


quadrática quando se exprime na forma V (x) = 
n 


= >. Jx X; X, onde q, = q, podendo-se 
nk=1 

representar por V (x) = x" Q x onde Q é a 

matriz simétrica dos coeficientes g,, . 

Para interpretar geometricamente a forma 
quadrática observa-se que no espaço RF. por 
exemplo, V (x) =x, +x4 corresponde a uma 
família de circunferências centradas na origem 
para diferentes valores V (x) = c'; do mesmo 
modo V (x) = 4, X1 + À, X5 com os parâmetros 
À, e à, positivos, representa um campo de elipses; 
e mais genericamente V (x) =q,, x; + 2q,,X, X,+ 
+ q Xs corresponde a uma família de elipses 


rodada em relação aos eixos coordenados, podendo 
2 


escrever-se V (X) = > T ix Xj X COM G5; = Go 
pk=1 

A forma quadrática definida positiva repre- 

senta o tipo de funções definidas positivas mais 

simples. Se os coeficientes q,. não são constantes 

mas dependem de x,, ..., x, a função 


n 
V' = > PR (PD 
pk=1 
define a pseudo-forma quadrática. 

A verificação da definibilidade de uma função 
efectua-se através do Teorema de Sylvester: A con- 
dição necessária e suficiente para que uma função 
V (x) na sua forma quadrática V (x) = x Q x 
seja definida positiva é que sejam positivos todos 
os sub-determinantes da diagonal principal da 
matriz Q 


n 


Ma) x, =x" 0 (x) x 


| 0 >0 =P = 


Jun, Ip he EE RC FT 
TJ, Jo | T 93 E dr Tn 


q qo 
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Também a forma pseudo-quadrática pode ser 
definida positiva se os coeficientes q, (X,,..-,X,) 
satisfazem ao teorema de Sylvester, isto é, se num 
ponto qualquer x, == O os valores dos coeficien- 
tes q, (X,) determinam a verdadeira forma qua- 


n 
drática > Ju (Xo) XX definida positiva, sendo 
k=] 
n 


portanto em X, o valor > Gu (X9) Xijo Xyo UM 
jk=1 
valor positivo. 

Para se verificar a definibilidade de uma dada 
função V (x) verifica-se se os determinantes 
D = qu. D, = q a D, = |Q| são 
Ta Too, 
todos positivos (função definida positiva). No caso 
em que — V (x) é definida positiva será V(x) defi- 
nida negativa, quer dizer, para provar que uma 
forma ou uma pseudo-forma V = x" Q x é defi- 
nida negativa utiliza-se o teorema de Sylvester 
para provar que — V == x" (— Q) x é definida 
positiva. 


3. DEFINIÇÕES DE ESTABILIDADE 


Considere-se um sistema autónomo [2] definido 
pela equação de estado x' = f (x) expressa nas 
variáveis de estado X. 

Por uma translação dos eixos coordenados do 
espaço de estado do sistema é sempre possível con- 
siderar um dado ponto de equilíbrio na origem 
do sistema, de entre os vários pontos de repouso 
que o sistema possa ter correspondentes a f(x) = 0. 

Na ocorrência de uma perturbação do sistema 
qué o retire transitóriamente do ponto de equi- 
líbrio, o sistema movimenta-se sobre uma trajtc. 
tória no espaço de estado, caracterizada pela 
variável tempo t. 

Admitida a origem O como sendo um ponto 
de equilíbrio do sistema, seja S (r) um domínio 
à volta da origem, de raio genérico r > o, cons- 
tituído por pontos satisfazendo a |X|| rr, que 
designamos por domínio inicial, no interior do qual 
se considera o início do movimento do sistema. 
O domínio S (r) pode ser arbitráriamente pequeno. 

De modo análogo se define o domínio de ingresso 
S (R), de raio genérico R > o, como sendo 
o maior dos domínios contendo a origem (ponto 
de equilíbrio) para o qual o movimento se veri- 
fica sempre no interior desde que tenha come- 
çado num dado domínio inicial com r<.R. É evi- 
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dente que o domínio de ingresso é função do 
domínio inicial, isto é, R é uma função de r. 

Nestas condições o ponto de equilíbrio (origem 
das coordenadas variáveis de estado) diz-se: 


a) Simplesmente estável (1) se cada trajectória ini- 
ciada num dado domínio inicial S (r) não 
sai para fora de um domínio de ingresso 
S (R), onde R ==R (r), quando t tende para 
infinito (fig. 1-a). 


e) 


Fig. 1 — Definição de estabilidade. 
a) Estabilidade simples 
b) Instabilidade. 
c) Instabilidade indiferente. 
d) Estabilidade assimptótica. 
e) Estabilidade assimptótica-global. 


b) Instável, em caso contrário (fig. 1-b). 

c) Indiferente, (?) quando existir uma trajectória 
em cujos pontos se verifica a estabilidade 
simples (fig. 1-c). 


(!) Também se designa na literatura por estabilidade em 


sentido de Liapunov. 


(*) Encontra-se na literatura a designação de estabili- 
dade fraca, 
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d) Assimptóticamente estável se cada movimento 
começado em S (r) não só permanece dentro 
de S(R) mas também se aproxima da ori- 
gem quando t tende para infinito (fig. 1-d). 

e) Assimptótica-globalmente estável se a estabilidade 
assimptótica corresponde à totalidade do 
espaço de estado (fig. 1-e). 


Note-se que a noção de estabilidade é caracte- 
rizada pela dimensão do domínio de ingresso, 
que no caso limite de ser infinitamente grande 
conduz ao conceito de estabilidade global. 


4. MOVIMENTO DE UM SISTEMA 


Para esclarecimento do conceito de estabili- 
dade analisa-se o movimento de um sistema no 
espaço de estado bi-dimensional, caracterizado 
por duas variáveis de estado x, e x, (o mesmo 
raciocínio generaliza-se para qualquer número de 
variáveis de estado, muma interpretação gráfica 
num hiper-espaço n-dimensional), quando se 
parte de um ponto inicial X, ligeiramente afas- 
tado da origem das coordenadas. Conforme a 
trajectória descrita pelo sistema, após este desvio 
transitório da origem admitida um ponto de 
equilíbrio, assim se terá um tipo de estabilidade. 

A fim de melhor se compreender as noções 
introduzidas compara-se cada um dos casos pos- 
síveis com uma evolução temporal relativa à pro- 
jecção do movimento no plano x,t de uma das 
variáveis de estado x. 


a) O sistema tende para um ponto estável: 


Caso 1 (fig. 2-a) — ponto estável final mais afas- 
tado da origem que o ponto 
inicial (estabilidade simples). 

Caso 2 (fig. 2-b) — ponto final mais perto da 
origem (estabilidade sim- 
ples). 

Caso 3 (fig. 2-c) — ponto final na origem (esta- 
bilidade assimptótica). 


b) O sistema tende para uma trajectória estável 
(estabilidade indiferente ou fraca): 


Caso 4 (fig. 2-d) — ponto inicial interior à tra- 
jectória estável. 

Caso 5 (fig. 2-e) — ponto inicial exterior à tra- 
jectória estável. 

Caso 6 (fig. 2-f) — ponto inicial scbre a tra- 
jectória estável. 
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c) O sistema não tem pontos estáveis (instabili- 


dade): 
Caso 7 (fig. 2-g) — ponto final no infinito. 


5. FUNÇÃO LIAPUNOV 


Como se indicou o método directo conduziu 
fisicamente ao estabelecimento de uma função 
«energia» que, satisfazendo à propriedade física 
de decréscimo continuo da energia quando t 
aumenta, permite apreciar a estabilidade. Estas 
características podem contudo exprimir-se de um 
modo genérico pela análise matemática, definin- 
do-se então a função Liapunov. 

Chama-se função Liapunov, em relação às variá- 
veis de estado X de um sistema dinâmico com a 
equação de estado x ==f(x), à função escalar 
V(x) definida positiva numa vizinhança da ori- 
gem e com derivadas parciais contínuas de modo 


o 9V 
a que a sua derivada V'(X) = >; 
dx 


i—l | 
V ()=f" grad V(x) seja definida negativa ou 
pelo menos semi-definida negativa nessa vi- 
zinhança. 

Com vista à interpretação geométrica deste 
conceito considere-se um sistema de 2.º ordem, 
definido por duas variáveis de estado, X = [x,, x,). 
Admita-se uma função qualquer V(x) =c, com 
c>0, onde o valor de c determina uma curva 
no espaço de estado designada linha de nível. 
Diferentes valores de c originam um campo de 
curvas de nível em todo o plano de estado (x,, x,). 
Supõe-se que a função V(x,, x,) é unívoca e é 
diferenciável em todos os pontos, existindo por- 
tanto em cada ponto sempre uma linha de nível. 
Observam-se contudo três pontos singulares 
(fig. 3): pontos extremos («pico» ou «fundo»), 
onde as linhas de nível se reduzem a pontos, 
e ponto «sela», no qual duas linhas de nível se 
contactam. Em qualquer destes casos o plano 
tangente à superfície V(x,, x,) é horizontal e por 


: dV dV 
conseguinte = 


dx, d X, 


f, (x) ou 


= 0, quer dizer, nestes 


pontos singulares será V' =0, 
Nestas circunstâncias, sendo V(x, x) uma 
função Liapunov, como se definiu, conclui-se: 


a) S6 pode haver pontos singulares na origem, 


Jd ) 
pois sendo V' = Ls E À f, definida 
dx, OR - 


2 
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Fig 2 — Tipos de movimento de um sistema: estabilidade simples (a, b), estabilidade assimptótica (c), 
estabilidade indiferente (d, e, f) e instabilidade (g). 


negativa será V' = 0Q na origem x=0 e no res- 
tante espaço onde é definida V' < 0. 

b) O ponto singular na origem só pode ser 
um ponto extremo visto que V'(x) deve ser 
continuamente monótonamente decrescente. 

c) A origem só pode ser um ponto fundo 
dado que a função V (x) é definida positiva, isto 
é, nula na origem O e positiva no restante espaço, 
pelo que a função Liapunov V(x) corresponde 
geometricamente a um «poço» (fig. 4- a). Assim 
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as curvas de nível à volta da origem são fechadas 
e os seus valores paramétricos crescem à medida 
que as curvas se afastam da origem. 

Ao espaço onde a função Liapunov é definida 
chamaremos domínio Liapunov. Dentro deste domi- 
nio uma trajectória deverá penetrar nas linhas 
de nível de fora para dentro (fig. 4-b) porquanto 
sendo V'< O ao longo da trajectória (só na ori- 
gem é V'==0) a função V(x(t)) terá de ser 
uninterruptamente decrescente. 
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D) 


Fig. 3 — Representação gráfica de uma função qual- 
quer V(x,, x9):1—pico, 2 - fundo, 3 — sela, 
a) Monte e poço da superfície V (x). 
b) Campo de linhas de nível. 


6. ESTABILIDADE ASSIMPTÓTICA 


Vê-se portanto que a função Liapunov, numa 
vizinhança da origem, corresponde a uma trajec- 
tória que se dirige para a origem. Nestas condições 
se existir essa função a estabilidade assimptótica 
fica assegurada. Observa-se todavia a hipótese 
admitida da vizinhança ser suficientemente 
pequena para que se verifique realmente a conver- 
gência da trajectória, isto é, embora o domínio 
Liapunov possa ser grande a região interessada 
limita-se ao chamado domínio de ingresso. 

Se o domínio Liapunov se estender a todo o 
espaço de estado, sendo V (x) —» «o quando a 
norma |'x|| —» «>, então também o domínio de 
ingresso corresponde à totalidade do espaço e 
caracteriza-se a estabilidade assimptótica global, 

Decisivo por conseguinte na apreciação da 
estabilidade assimptótica é encontrar uma função 
Liapunov que determine o maior domínio de 
ingresso possível. Deste modo se reparte a noção 
de estabilidade assimptótica em pequeno, em 
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a) 


b) 


Fig. 4 — Representação gráfica de uma função Liapunov 
Vixp X9): O — ponto fundo (origem). 
a) Poço da superfície V (x). 
b) Campo de linhas de nível e trajectória 
numa vizinhança da origem. 


grande e em total, conforme as dimensões do 
domínio de ingresso. 

A estabilidade assimptótica em pequeno é caracteri- 
zada por um domínio de ingresso de reduzidas 
dimensões, podendo portanto praticamente cor- 
responder a uma instabilidade se esse domínio 
for excessivamente reduzido. Por exemplo, um 
sistema em equilíbrio a 100 V se possuir um 
domínio de ingresso com o raio de Im V é 
praticamente instável dada a possibilidade normal 
de flutuações da tensão acima de 0,001º/,. 

Na estabilidade assimptótica em grande define-se 
um certo domínio de ingresso limitado no espaço 
de estado (contendo a origem e onde a função 
Liapunov é definida, verificando-se por conse- 
quência em qualquer ponto desse domínio a 
estabilidade assimptótica), podendo ser apenas 
parte do domínio Liapunov, pelo que na restante 
zona do domínio Liapunov se poderá verificar a 
instabilidade pois a trajectória, dirigindo-se dentro 
do domínio Liapunov, onde V' < 0, sempre para 
as linhas de nível de menor parâmetro, poderá 
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eventualmente atingir a zona do espaço em 
que V/> 0 (exterior ao domínio Liapunov) sem 
se dirigir por consequência para a origem, isto é, 
sem penetrar no domínio de ingresso (fig. 5). 

No caso de estabilidade assimptótica em total o 
domínio de ingresso pertence à totalidade do 
espaço de estado, correspondendo portanto à 
estabilidade assimptótica global, para o que se 
exige que a função Liapunov seja infinita quando 
o ponto X do espaço de estado se encontra no 
infinito: V (x) — co para ||X || —» os, 


so: 
Mo 

/ f Ei | ai 
[Fo | 
0107/44 
ns da 
o 
| 

V>0 | v<0  VSO 


Fig. 5 — Estabilidade assimptótica em grande (domínio 
Liapunov para V' < 0). 


7. ESTABILIDADE INDIFERENTE E INSTA- 
BILIDADE 


Uma função Liapunov é definida por duas 
condições na vizinhança da origem: 12) V (x) 
definida positiva; 22) V' (x) definida negativa. 
Então o ponto de equilibrio é assimptóticamente 
estável. 

O facto de V (x) ser definida positiva determina 
na origem um ponto fundo e no espaço de estado 
linhas de nível V (x)=-c, com c positivo, de 
parâmetro c crescente em todas as direcções a 
partir da origem (envolvendo este ponto). 

Sendo V” (x) definida negativa, V' (x) “0 em 
todos os pontos x + O da vizinhança da origem 
e portanto V(x (t)) decrescente sucessivamente 
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no tempo, quer dizer, a trajectória x (t) tende 
para a origem onde V' (x) == 0. Este processo 
caracteriza a estabilidade assimptótica. 

Mas se V' (x) é semi-definida negativa, então 
Vº (x) <' O na vizinhança da origem, havendo 
portanto pontos, além da origem, onde V' (x) =0. 
Como V(x (t)) não pode crescer a trajectória x (t) 
não se pode dirigir para as curvas de nível de 
maior parâmetro permanecendo então sobre a 
linha de nível correspondente ao início do movi- 
mento ou noutra inferior. Neste caso a origem 
é estável dizendo-se que se está em presença de 
estabilidade fraca (fig. 6-b), que preferimos referir 
por estabilidade indiferente, pois deslocando o sistema 
para qualquer dos pontos dessa linha de nível o 
sistema pérmanece em equilíbrio. 

É evidente que se V' (x) for definida positiva 
a origem, para a correspondente função V (x), é 
instável. De facto as curvas de nível continuam 
a circundar a origem com parâmetros crescentes 
em todas as direcções mas V' (x) > 0 exterior- 
mente a O, quer dizer, V (x (t)) cresce no tempo, 
de modo que a trajectória X (t) se dirige perma- 
nentemente para as linhas de nível com maior 
parâmetro, o que conduz à instabilidade da origem. 

Todavia o facto de ser V' (x) definida positiva 
significa também que a função arbitrada V (x) 
não é uma função Liapunov. Por outro lado se 
a origem não for assimptóticamente estável não 
existe nenhuma função Liapunov. Conclui-se 
assim que a existência de uma função Liapunov 
é condição suficiente mas não necessária de 
estabilidade assimptótica. Por outras palavras: 
se não se encontrar uma função Liapunov não 
se pode afirmar a instabilidade. 


8. CRITÉRIO DE ESTABILIDADE 


A discussão efectuada permite estabelecer um 
critério de verificação da estabilidade, que designa- 
remos por crilério Liapunoo por se basear na 
formação de uma função Liapunov. 

Num sistema autónomo de equação de estado 
x == f (x) com um ponto de equilíbrio na origem 

= 0, se existir uma função V(x) definida 
positiva e diferenciável muma vizinhança da 
origem onde 


a) V'(x) é definida negativa, o ponto de 


equilíbrio é assimptoticamente estável (fig. 6-a) 
e conforme a dimensão do domínio de 
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ingresso assim será assimptoticamente está- 
vel em pequeno, em grande ou em total 
(global). 

b) V' (x) é semi-definida negativa, o ponto de 
equilíbrio é fracamente estável ou de estabilidade 
indiferente (fig. 6-b). 

o) V' (x) é definida positiva, o ponto de equi- 
líbrio pode ser instável (fig. 6-c). 


Este critério é facilmente estendido a sistemas 
de tempo não invariante ou não estacionários 
com a equação de estado x” = f (x, t). 


Fig 6 — Informação dada por uma função definida 
positiva V (x). 
a) Estabilidade assimptótica (V' definida ne- 


gativa). 

b) Estabilidade indiferente (V' semi-nega- 
tiva), 

c) Possibilidade de instabilidade (V' definida 
positiva). 


374 


9. ESTABILIDADE ABSOLUTA NUM SEC- 
TOR 


Considere-se um sistema de controle com um 
elemento não linear de característica F (7), exci- 
tado pela grandeza 7, e um elemento linear de 
característica G (7) em relação à variável 7 (fig. 7-a). 

As condições impostas pelo teorema Sylvester 
na definibilidade de V e V' definem um sector de 
estabilidade [k,, k,] caracterizado genericamente 
por k, <. F(9)/2 <. k, contido no 1.º e no 3.º 
quadrantes do plano (F (7), 7) devendo todas 
as características do sistema estar contidas nesse 
sector e passando pela origem para que o ponto 
de equilíbrio seja assimptóticamente estável 
(fig. 7-b). Nestas condições define-se a estabilidade 
absoluta num sector. 

Como é evidente o ângulo do sector depende 
da correspondente função V formada pelo que 
no caso de diferentes métodos de formação da 
função Liapunov conduzirem a funções V dis- 
tintas também os respectivos sectores serão em 
geral diversos. Dado que a existência de uma 
função Liapunov é uma condição suficiente mas 
não necessária de estabilidade conclui-se que 
é válida a lei da sobreposição dos sectores, pois 
cada um corresponde a uma parte da zona com- 
pleta de estabilidade. 


Hig.7 — Estabilidade absoluta num sector de um sis- 
tema de controle. 
a) Esquema de blocos. 
b) Sector de estabilidade, 
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10. EXEMPLO DE PROVA DA ESTABILIDADE 


Para exemplificação dos conceitos discutidos 
considere-se um sistema constituído por três ele- 
mentos: uma mola de constante k,, de maneira 
a que a força exercida para um desvio x é, 


F, = — k, x onde o sinal — significa fisicamente 


Fig. 8 


oposição ao desvio, um amortecedor de constante 
k,, sendo anaálogamente F, = — k, x” e uma 
massa m (fig. 8-a). 


10.1. Pelo método directo de Liapunov 


Aplique-se o método directo, de Liapunov, 
a este sistema. 
A equação do movimento do sistema será 


e di a! o dao DR 


donde, pela introdução das variáveis de estado 


se obtém a equação de estado 
, k, ko 


16d e o E e DO E 
2 m I m 2 


O ponto de equilíbrio verifica-se portanto 
para x, == 0 ou seja para x, == x, = 0. 
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A energia do sistema é igual à soma da energia 
cinética com a energia potencial, sendo 


1 ps 
E = mx = mx; 
o o 
e E = | F,dx = | —yxde= 1x7 
x X1 


donde E == 


função que é positiva em todo o espaço de estado 
R? e nula na origem (por outras palavras, E é 
uma função definida positiva). 

Quando o sistema sofre um desvio da sua 
posição de equilíbrio a energia varia no tempo de 
acordo com a derivada 


+ f ' 
ou seja E =k,x,x,+t+mx,x, 


e, pela definição de x, e com a equação de estado, 


k, k, ) 
L = kpxoxa timão ( — Lx a 
obtendo-se E =—k,x) 


quer dizer, a energia diminui constantemente em 
todo o espaço de estado tendendo o sistema 
necessáriamente para o repouso (note-se que E” 
é uma função definida negativa). 

Conclui-se assim, sem resolver a equação de 
estado, que a origem é assimptótica-globalmente 
estável. 

Mas para sistemas mais complexos o método 
directo apresenta duas dificuldades: é difícil 
exprimir a energia e, se se consegue obter uma 
expressão, é difícil analisar a variação temporal 
da energia (muitas vezes não decresce monótona- 
mente). 


10.2. Por uma função Liapunov (critério 
Liapunov) 


O critério Liapunov de investigação da estabi- 


lidade consiste na determinação de uma função 
Liapunov. 
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O método directo sugere as variáveis de estado 
X, Et =; com Xs =-—a,Xxj—açx,e a fun- 
ção V (x) = x] + x) que deverá satisfazer às con- 
dições Liapunov. Mas, apesar da função V(x) = 
4) + K ser definida positiva, a sua derivada 
Vetj==?(a—Ie xa ax não é defi- 
nida negativa nem semi-definida negativa. Resulta 
que a função arbitrada não prova a estabilidade 
do sistema. De facto as curvas de nivel seriam 
circunferências centradas na origem e a trajectó- 
ria intersectaria as linhas de nível de dentro para 
fora (fig. 8-b). 

Admita-se porém VM) =q,X, +H2G,X,X, + 
Moss de coeficientes indeterminados, definindo 
uma família de elipses como linhas de nível, 
função que deve ser definida positiva. Calculem- 
-se de acordo com o teorema Sylvester, os coefi- 
cientes q,,,9],, €Gy demodo aque V* (x) — e: — 
Rj , portanto definida negativa. Os valores dos 
coeficientes assim calculados verificam realmente 
as condições impostas pelo teorema Sylvester à 
função V(x) e portanto V(x) constituiu uma fun- 
ção Liapunov, concluindo-se a estabilidade assimp- 
tótica do sistema, 

Na verdade é fácil verificar pela derivação de 
Ví(x) e identificação com V'(x) =—*) — Xo que 
deverá ser 


1 
qo — 
2a, 
l+a, 
É Do 
a 2a,a, 
ão + a, (1 + a,) 
O a id 


2a,a, 


valores que verificam as condições dadas pelo 
teorema Sylvester à função V(x) 


%120 
9 
1% — 9220 


pois por substituição de valores em q,, Go — 


— GD se deduz que terá de ser a, >0 o que 
conduz à verificação de q, 0 quando a, >. 


11. FORMAÇÃO DE FUNÇÕES LIAPUNOV 


O problema da estabilidade segundo o critério 
Liapunov depende da obtenção de uma função 
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Liapunov, nomeadamente de uma «boa» fun- 
ção Liapunov, que determine o maior domínio 
de ingresso possível. 

Não existe porém nenhum método geral que 
conduza à sua formação no caso de sistemas não 
lineares. Nos casos práticos dos sistemas de 
controle existem vários métodos de entre os quais 
se apresentam algumas ideias basilares. 


11.1. Método em sistemas lineares 


Num sistema linear, traduzido pela equação 
de dinâmica x'==Ax, onde À é uma matriz 
constante (independente de X), considere-se uma 
forma quadrática V(X) — x 0x positiva para 
x=0 e V(0) =-0, onde Q é uma matriz simé- 
trica real definida positiva, sendo portanto 
Vogy=x "ax + x" Qx. Como x =Ax exT = 
=x A! será V'(xg=x' (AQ + A)x, donde 
V'(g=—x"PxcomA!QLQOA=-—P. Deste 
modo se P for definida positiva será V'(x) defi- 
nida negativa, quer dizer, o sistema é estável. 

O problema resume-se pois a calcular os coe- 
ficientes q, da matriz Q através de ATO + 
+ GA = — P, tendo-se feito P definida positiva, 
o que permite formar uma função Liapunov 
através de V (x) =x" Qx. 

Com este raciocínio o critério de estabilidade em 
sistemas lineares pode enunciar-se: É condição 
necessária e suficiente para que um sistema linear 
de equação dinâmica x ==A x seja estável que 
para uma dada matriz simétrica definida positiva 
(real) P exista uma matriz simétrica definida Q 
tal que A Q-+-0A=—P, 

Como V(x) resulta uma função Liapunov é 
evidente que a condição é suficiente. Pode-se 
também demonstrar que é condição necessária. 


11.2. Método aproximado em sistemas linea- 
rizados 


No caso de sistemas não lineares a matriz À 
da equação dinâmica x =-A(x)x deixa de ser 
constante, dependendo pelo menos um dos seus 
elementos a,, do vector de estado x. 

Esta matriz À (x) pode contudo decompor-se 
em duas sub-matrizes, A, constante e A, (x) 
dependente de x, de modo que A (x) = A, + A (x), 
sendo portanto x = À, x + Ap (x) x que define 
a linearização em torno da origem se A, (x) O 
quando x — O. Esta condição não é indispensável 
para se definir um método aproximado. 
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Partindo da forma quadrática V (x) = x 0x, 
com os coeficientes indeterminados, será V (x) = 
=x PxcomP(xg)=a" Q+0A que se pode 
decompor em duas partes, uma matriz P, cons- 
tante e outra Pp, (x) dependente de x, tal que 
P(x) =P, + P, (x), sendo por conseguinte 
PL=A,/0+0A eP(x)=Aç(X)0 +04, (x). 
Fazendo P, = — |, onde | é a matriz unidade; 
pode calcular-se Q através de A, Q+04,=—| 
pois esta equação é resolúvel em Q desde que a 
matriz À, tenha apenas valores próprios (raizes 
do determinante det (A, — 4 |=0) com parte 
real negativa, isto é, sempre que a matriz À, 
corresponda a um sistema linear estável. Nestas 
condições a forma quadrática V = x" Q x é defi- 
nida positiva. 

Por outrolado P(x) =— | + P, (x) ou P (x) = 
=— | + AR (x) 0 + Q Ap (x) corresponde 
à forma quadrática V' — x" P (X) x pelo que por 
meio do teorema de Sylvester se pode calcular o 
domínio onde V” é definida negativa. 

Note-se que na medida em que A, (x) 0 
quando x» O também Pp, (x) - O e portanto na 
vizinhança da origem 0 é P (x) =—l isto é, 
V'x)=—x" Ixou seja ainda V' (x) =x) Eu + 
donde se vê que V” é definida negativa na vizi- 
nhança da origem 0, 


11.3. Método de Aiserman 


Baseado nas considerações anteriormente ex pos- 
tas Aiserman [3] substituiu a característica do 
elemento não linear do sistema em estudo por 
uma característica linear aproximada, o que pelo 
critério da estabilidade dos sistemas lineares 
permite determinar os coeficientes q,, da matriz Q 
e portanto definir uma função V que se toma 
como função Liapunov do sistema pela introdução 
da característica não linear. 


11.4. Método de Szego 


O modo como Szego [4] restringiu a função V' 
conduziu a supor que na função V os coeficientes 
qi; são funções apenas de x, e x,, quer dizer, 
Giu (Xj, X,) e não g,;, (X). Os coeficientes a,. sendo 
considerados polinómios em x, e x,, não são 
função de x,, pelo que V” é sempre uma equação 
algébrica do 2.º grau em x, e portanto a equação 
V'==0 pode ser resolvida em x, pela forma 
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quadrática. A resolução desta equação define 
duas superfícies no espaço de estado mudando 
o sinal de V” quando estas superfícies se cruzam, 
quer dizer, entre as duas superfícies os pontos 
do espaço são caracterizados por exemplo por 
Vº<L 0, e nos pontos da zona exterior às super- 
fícies por V' 0, sendo evidentemente sobre 
elas V' = O (fig. 9). Se se força estas duas super- 


Fig. 9 — Método de Szego no espaço R2, 


fícies a serem coincidentes, V' não mudará de 
sinal em todo o espaço de estado e estabelecem-se 
as condições de Liapunov. 


11.5. Método de Schultz e Gibson 


No método de Aiserman partindo-se da forma 
quadrática V com coeficientes indeterminados 
calcula-se o vector grad V, que permite formar 
V'=f grad V e determinar os coeficientes de 
maneira a que V seja definida positiva. O método 
de Schultz e Gibson, designado método de 
gradiante variável [5], procede inversamente : 
parte do vector grad V, cujas componentes são 
formalmente lineares de coeficientes indetermi- 
nados, escolhendo-se os coeficientes q, de tal 
modo que V seja definida positiva e V' definida 
negativa ou semi-definida negativa. 

Em geral este método é mais fácil de aplicar 
nos casos práticos do controle automático, depen- 
dendo contudo da escolha dos coeficientes q, 
que no método de Aiserman se fixam pelo 
desenvolvimento do cálculo. 

Assim se no sistema de equação x =x (f) 
com f ()) = O existe, numa vizinhança da origem, 
um vector real grad V diferente de zero excepto 
d ( OV ) no d o) 


dx; rr OX NOx, 


na origem tal que 
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e se V' for definida 


x 
” ) 
- | Es dx sos + Edi 
' X1 n 


negativa com V = 


dx ) definida 
o d x n 


o 
/ 


positiva, o sistema é assimptôticamente estável 
nessa vizinhança da origem. 

O método do gradiente variável procede 
similarmente ao método de Szego dando contudo 
maior âmbito à função V através da introdução 
de integração na sua determinação. 


11.6. Metodo de Lurie 


Lurie e Letov admitiram um integral na definição 
de V, «forma quadrática mais um integral sobre 
a não linearidade», incluindo assim em V um 
maior número de termos, o que dá a q, maior 
generalidade. A dificuldade resultante deste maior 
âmbito consiste nas restrições de V” para que 
seja pelo menos semi-definida negativa. 

Decompondo o elemento linear [6], equiva- 
lente a todos os elementos lineares do sistema, 
em partes caracterizadas por c, e 2, (fig. 10), que 
definem respectivamente o vector C e a matriz A 
com os valores próprios 4, = «,, a equação diná- 
mica do elemento linear x == À X será represen- 


Fig. 10 


tada por Xº == AX e portanto, considerando 
a característica f (2) do elemento não linear, 
a equação dinâmica do sistema resulta x = À x + 
-+- e f (5), sendo € um vector unitário. 


Pelo diagrama de blocos 7 = — y com y = c! x 
ou seja 7 = — c! x donde se deduz 
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o =bTx-rf(:) comb! = - c'A 


a 
e [= €= 2,€C. 


Lurie parte da função 


d - TG, IG Xi; Xx : 
V=x— ——— t(g) dE, 
2 2 A+ ds +] ida 


=| k= 


que é definida positiva quando o integral for posi- 
tivo, isto é, para f(2)/2 > :7>0, ou seja, se a 
característica do sistema não linear está no 1.º e 3.º 
quadrantes do plano (f(7),7). Para que a partelinear 
do sistema seja estável admite-se que a parte real 
de », é negativa. Nestas circunstâncias deduz-se 


1 2 n 
V=—rf (o)—( ar] + f (7) Dm 
vj º i=1 


/ n 
so % 
A b;— 29; é, ET que será definida nega- 
k=1 0 "4 *k 
n 

a b > o ceu 
tiva quando r > 0 e b= 27], 1), pois 
se anula a única parcela positiva de V”. Satisfei- 
tas estas condições para i==1,2,...,n o sistema 


é assimptoticamente estável. 
11.7. Método de Krasovski 


Seja um sistema x'=f ix) com o ponto de 
equilíbrio x = O e em que f (x) admite derivadas 
parciais em ordem a x, (comi=1,2,..., n). 

o fo 
Considerando a matriz Jacobiana J (x) = Ê 1 
dx 
k - 


+ 


ce a matriz simétrica ] (x: for definida negativa 
(condição suficiente) o ponto de equilíbrio x = O 
éestávele V (x) = f" (x). f (x: é uma função 
Liapunov. Se além disso E” (x) . f(x) —co quando 
| X| — co então o ponto de equilíbrio é assimp- 
tótica-globalmente estável [7]. 


12. CONCLUSÕES 


Tendo-se apresentado um exame exaustivo dos 
tipos de estabilidade possíveis para um sistema 
definiu-se o critério de estabilidade pela definição 
de uma função Liapunov. O problema reduziu-se à 
determinação desta função, tendo-se discutido 
vários métodos de formação. 
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4 HERÓIS DA LUBRIFICAÇÃO > 


Merelim 


Há uns 1400 anos, mais coisa menos coisa, um sujeito 
que respondia ao nome de Artur arrancou uma espada de 
uma pedra e, assim, sem mais nem menos, tornou-se Rei de 
Inglaterra. Mas quem lhe facilitou a tarefa, não foi outro 
senão um mágico da córte que se chamava Merelim. 

No século VI, a Inglaterra lutava com uma grande 
falta de reis. Por isso, alguém se lembrou de criar um teste 
de perícia que consistia na aparentemente estranha opera- 
ção de retirar uma espada de uma bigorna, assente sobre 
uma grande pedra. O anúncio publicado rezava assim: 
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Merelim 


«Quem conseguir arrancar esta espada 
será proclamado Rei de Inglaterra!» 


As tentativas que diáriamente faziam os nobres can- 
didatos ao trono vago redundavam em fracasso. 

Artur, o nosso homem, que fazia uns biscates por 
aqui e por ali, ansioso que estava — coitado! — de arran- 
jar um lugar com mais futuro e maiores possibilidades, con- 
sultou o seu amigo Merelim, mágico da córte. 

Merelim estudou a situação e meteu mãos à obra. 
Tudo ele tentou para conseguir resultados positivos: os tri- 
viais passes de magia, os pózinhos de perlimplim e, até, 
algumas maldições bastante fortes, diga-se. Mas a espada 
não se mexia. 

Ao cabo de tantas e tão infrutiferas tentativas, Mere- 
lim disse: «O problema consiste na lubrificação e a minha 
magia não resulta. Só com ensaios de laboratório podere- 
mos chegar a uma conclusão qualquer». 

«Mas a minha prova é daqui a dois dias!», protestou 
aflito, o nosso Artur. «Não temos tempo para isso». 

Merelim, contudo, informou-o de que o seu fornece- 
dor de lubrificantes tinha laboratórios espalhados por todo 
o mundo e que a resposta desejada poderia ser obtida 
muito depressa. 

Se assim o pensou, melhor o fez. Sem perda de um 
minuto, o problema foi entregue ao laboratório da MOBIL, 
que tinha instalações no próprio bairro. 

Os técnicos iniciaram imediatamente os ensaios que 
se impunham. Ao cair da noite, tinham chegado a uma 
conclusão positiva: uma mistura especial de óleo de asas 
de morcego e de gordura de escorpião, mistura esta que 
faria deslizar a mais pesada das espadas que estivesse enta- 
lada na mais teimosa das pedras. 

Graças à rapidez dos serviços laboratoriais da MOBIL 
e à oportunidade com que Merelim, um mago que soube 
reduzir-se à sua modéstia, desistindo da magia quando ela 
não servia, aplicou a massa lubrificante «TESE» (Tira Espa- 
das Sem Espinhas) da MOBIL, Artur conseguiu o emprego 
que pretendia: reinar. 

Os clientes de hoje também confiam na MOBIL para 
os lubrificantes e serviços de que necessitam. 

A MOBIL dispoe de laboratórios espalhados por todo 
o mundo. Isto significa que há sempre um laboratório perto 
de qualquer problema de lubrificação, um laboratório equi- 
pado para o resolver pelo método mais rápido e seguro. 

Este é mais um dos serviços específicos que se tra- 
duzem por lubrificações excepcionais à disposição dos com- 
pradores dos nossos produtos. 

Nao é magia. E MOBIL! 
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PRDOREDIOM - 


a) O método directo de Liapunov prova à 
estabilidade sem se resolver a equação diferencial 
de estado de sistema, sendo encarado como um 
caso particular do critério de estabilidade. 

b) O critério Liapunov permite avaliar a esta- 
bilidade de um modo não unívoco, pois em prin- 
cípio pode construir-se uma infinidade de fun- 
ções Liapunov (as condições de estabilidade do 
critério são apenas suficientes e não necessárias). 

c) O domínio de ingresso na estabilidade em 
grande depende da função V (x), sendo difícil a 
formação da função Liapunov mais adequada. 

d) Para sistemas até 3.º ordem é possivelumainter- 
pretação geométrica, e para sistemas de ordem supe- 
rior é muito difícil a formação de funções Liapunov. 

e) Existem diversos métodos para formar fun- 
ções Liapunov mas nenhum ideal. 

Na prática dos sistemas de controle o critério 
Liapunov perdeu parte do seu interesse, e portanto 
todos os métodos discutidos de formação de fun- 
ções Liapunov com o critério de Popov, baseado 
na função de frequência modificada. 
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